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We use the continuous-time interaction expansion (CT-INT) quantum Monte Carlo method to 
calculate the phonon spectral function of the one-dimensional Holstein-Hubbard model at half¬ 
filling. Our results are consistent with a soft-mode Peierls transition in the adiabatic regime, and the 
existence of a central peak related to long-range order in the Peierls phase. We explain a previously 
observed feature at small momenta in terms of a hybridization of charge and phonon excitations. 
Tuning the system from a Peierls to a metallic phase with a nonzero Hubbard interaction suppresses 
the central peak, but a significant renormalization of the phonon dispersion remains. In contrast, the 
dispersion is only weakly modified in the Mott phase. We discuss finite-size effects, the relation to 
the dynamic charge structure factor, as well as additional sum rules and their implications. Finally, 
we reveal the existence of a discrete symmetry in a continuum field theory of the Holstein model, 
which is spontaneously broken in the Peierls phase. 

PACS numbers: 71.10.Pm, 71.45.Lr, 71.30.+h, 71.38.-k 


I. INTRODUCTION 

Electron-phonon interaction plays a crucial role for the 
physics of materials [1], and gives rise to phenomena such 
as the Peierls instability [2], superconductivity, and po- 
laron formation [3]. Recently, interest in electron-phonon 
interaction has also been boosted by experiments involv¬ 
ing photo-induced phase transitions between insulating 
Peierls and metallic states [4], see Ref. [5] for a review. 

The study of microscopic electron-phonon models has 
a long and rich history. Much of the recent progress on 
the theoretical side resulted from the development of ex¬ 
act numerical methods, most notably the density-matrix 
renormalization group [6-8], and quantum Monte Carlo 
(QMC) methods [9-15]. Whereas the case of a single 
electron (the polaron problem) can be solved to machine 
precision using a variational basis construction [16], fi¬ 
nite band fillings are still a challenge. Here, we discuss 
one-dimensional (ID) models. 

Of particular interest with regard to experiment is the 
calculation of spectral properties such as the electronic 
spectral function that may be compared to angular- 
resolved photoemission spectra. For the spinless Hol¬ 
stein model at half-filling, it reveals the opening of a 
Peierls gap for strong electron-phonon coupling [17-19], 
the experimentally observed shadow bands [20] arising 
from the new periodicity of the lattice [21], and soli- 
ton excitations [21]. The phonon spectral function and 
the dynamic charge structure factor—experimentally ac¬ 
cessible via neutron scattering—reveal, for example, the 
softening of phonon excitations near the Peierls transi¬ 
tion [19, 22, 23], and distinguish soft-mode behavior in 
the adiabatic regime from central-peak behavior in the 
nonadiabatic regime. For a single electron, the phonon 
spectrum can again be calculated to arbitrary precision 
[24]. The phonon spectral function and the renormal¬ 
ized phonon frequency of the spinless Holstein model 
were obtained with the projector-based renormalization 
method [17, 23, 25]. The electronic spectral function of 


the Holstein-Hubbard model was calculated to character¬ 
ize the metallic, Mott, and Peierls phases [26-28], and to 
study the impact of electron-phonon coupling on spin- 
charge separation [29, 30]. A more complete review of 
previous work can be found in Refs. [21, 28]. 

Here, we show that the CT-INT method, which was 
applied to a number of electron-phonon problems [15, 
21, 28, 31-34], can also be used to calculate the phonon 
Green function, and thereby the phonon spectral func¬ 
tion. The method is free of a Trotter error, and does not 
require a cutoff for the phonon Hilbert space. We use 
it to calculate the phonon spectral function of the half- 
filled Holstein-Hubbard model. The results are discussed 
in the context of previous work on these and other mod¬ 
els, and additional insights are provided with the help of 
the random phase approximation and field theory. 

The paper is organized as follows. In Sec. II we define 
the model. In Sec. Ill we discuss the method. Our re¬ 
sults are presented in Sec. IV, followed by a discussion in 
Sec. V. Section VI contains our conclusions, and the Ap¬ 
pendix discusses the relation between charge and phonon 
spectra, sum rules, and their implications. 


II. MODEL 

We consider the Holstein-Hubbard model 

H=-tJ2 (4A+1.T + H - c -) +U J2 (!) 
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The first term describes the hopping of electrons be¬ 
tween neighboring lattice sites with amplitude t. The 
second term captures the repulsion between electrons at 
the same lattice site; for U = 0, Hamiltonian (1) becomes 
the spinful Holstein model [35]. The lattice is described 
in the harmonic approximation, with the displacement 
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(momentum) at site i given by Qi (Pi). The bare op¬ 
tical phonon frequency is uoo = \/K/M. The electron- 
phonon interaction is of the density-displacement type, 
with pi = hi - 1, hi = ^2 a h ia , and h i<7 = c\ a c ia . The 
dimensionless ratio A = g 2 / (AKt) is a useful measure of 
the electron-phonon coupling strength. 

We studied Eq. (1) at half-filling, corresponding to 
(hi) = 1, on chains with L sites and with periodic bound¬ 
ary conditions. We use t as the unit of energy, and set 
the lattice constant, M, and h to one. 

The Holstein-Hubbard model captures the Peierls tran¬ 
sition driven by strong electron-phonon coupling, the 
Mott physics related to strong electron-electron interac¬ 
tion, and the competition between these two phenomena 
when U ~ A [26, 36-39]. Although some open ques¬ 
tions remain regarding the phase diagram and the low- 
energy physics, the prevailing picture is that the model 
supports an intermediate metallic phase with a spin gap 
but gapless density fluctuations [28, 39], which is adi- 
abatically connected to the metallic phase of the spin¬ 
ful Holstein model (U = 0) [8]. The metallic behavior 
for A < \ c (U,u)o) is the result of quantum lattice fluc¬ 
tuations that destroy the long-range Peierls order. Its 
extent therefore depends on the phonon frequency: For 
uoo = 0 a Peierls state exists as soon as 4At > I/, whereas 
for cj 0 = oo a Peierls state is completely absent from the 
phase diagram [36-39]. In the classical limit and for half- 
filling, the charge order in the Peierls phase corresponds 
to spin-singlet pairs of electrons (singlet bipolarons) re¬ 
siding on every other lattice site. 

We are not aware of any previous results for the phonon 
spectral function of the ID Holstein-Hubbard model. 
The phonon spectrum of the ID, half-filled spinless Hol¬ 
stein model was studied in Refs. [17, 19, 22, 23, 25]. Nu¬ 
merical results were obtained using exact diagonalization 
(restricted to small cluster sizes) [19], cluster perturba¬ 
tion theory (with artifacts related to the inherent trans¬ 
lation symmetry breaking) [19], and QMC (restricted 
in cluster size and temperature, and with a Trotter er¬ 
ror) [22]. Regarding analytical work, the projector-based 
renormalization group approach [17, 23, 25] provides the 
most comprehensive results. For ujq < £, the renormal¬ 
ization of the phonon spectrum at large q was studied 
by means of the dynamic charge structure factor [21, 28]. 
A general feature of Holstein models is that for A = 0, 
the phonon dispersion uj q = cj 0 is trivial. However, 
for nonzero couplings, a renormalization takes place and 
ujq i —y ujq . For the spinless Holstein model, uj q has been 
calculated in Refs. [17, 23]. 

III. METHOD 

The CT-INT QMC method permits one to simulate 
rather general fermionic actions (including retarded and 
nonlocal interactions) [40]. Exploiting the fact that the 
phonon degrees of freedom can be integrated out exactly 
to obtain a fermionic action [15, 41], it was applied to a 


number of different electron-phonon problems [15, 21, 28, 
31-34]. The method relies on the imaginary-time path- 
integral formulation of the partition function, which is 
calculated using a weak-coupling perturbation expansion 
[15]. A strong-coupling formulation also exists, but is less 
suited for lattice problems [42]. 

For the model (1), the path-integral representation of 
the partition function takes the form 

Z = (V{c,c) e-SoM-SiM [ V(q) e ~ s ^’ c ^ > ( 2 ) 


where we used the coherent-state representation Ci a \c) = 
Ci a \c) with Grassmann variables Ci a for the fermions, 
and the real-space representation Qi\q) = qi\q) for the 
phonons. We split the action into the free-fermion part 
So , the Hubbard interaction Si, and the remainder S ep 
containing the free-phonon part as well as the coupling 
of the displacement fields to the electrons, given by 

SeV = !o rfr F{^T 2 (T) + - 9 q i(T)Pi(T)} ■ (3) 


Integration over the fields q in Eq. (2) leads to an effective 
fermionic action S 2 that can be simulated with the CT- 
INT method in the same way as fermionic models [15]. 

Previous applications of the CT-INT method were re¬ 
stricted to correlation functions of fermionic operators. 
By using a generating functional with source fields that 
couple to the lattice displacement fields, we can derive an 
estimator for the phonon propagator D^(r) in terms of 
the time-displaced charge correlation function. Explic¬ 
itly, to measure the phonon propagator 

Ai(r) = (ft(r)^(O)), (4) 

we add a source term S v = — f dr JT Vi( T )li( T ) to Sep, 
where r]i(r) is a real field. After integrating out the dis¬ 
placement fields, we arrive at an effective action describ¬ 
ing the electron-phonon interaction, 

S * = -^2JJ 0 drdT 'Yl IN r ) +9~ 1 Vi(' r )\ (5) 

x D°(t - t') [Piir') + g^rnir’)} . 
The appearance of the free phonon propagator 



defined for —/3 < r < yd, leads to a retarded interaction 
in imaginary time. The interacting propagator D^(r) 
can be obtained from Eq. (5) via a functional derivative 
with respect to gpr) in the limit g —» 0. The result, 


D ij (T) = D°(T)6 ij +g 2 Jf dTidT 2 D°(r 1 - t) (7) 

x D°(t 2 ) (Pi(ri)ft(r 2 )) , 
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FIG. 1. (Color online) Phonon spectral function B(q,uj) of the spinful Holstein model. Dashed lines correspond to ujo — 0.5 1. 
Here, U — 0, and f3t = L = 42. The Peierls transition occurs at A c ~ 0.25 [37, 38]. Color scheme based on Ref. [43]. 


is the sum of the free propagator and an interaction term 
that involves the time-displaced charge correlation func¬ 
tion via a double convolution with D°(r). The charge 
correlation function can be accurately measured with the 
CT-INT method. Equation (7) has previously been de¬ 
rived for the Anderson-Holstein impurity model [44]. 

We evaluated Eq. (7) for each bin average of 
(ni(r)rij( 0 )), and used a binning analysis to obtain re¬ 
liable statistical errors. To minimize systematic errors 
in the numerical integration, we carried out one of the 
integrals in Eq. (7) analytically by exploiting the peri¬ 
odicity of D°(r) and (ni(r)rij( 0)). Measurements were 
made on an equidistant grid on the imaginary time axis 
with Art = 0.1. From Dij(r) we obtained the phonon 
spectral function (we define A nm = E n — E m ) 

B(q,u) = \E \{ m \Qq\ n )\ 2e ~ l3EmS (^ ~ A nm) ( 8 ) 
n,m 

by carrying out a Fourier transformation and using the 
stochastic maximum entropy method [45] for the ana¬ 
lytic continuation. The use of pi(r) instead of n^(r) in 
Eq. (7) amounts to subtracting the T = 0 static contri¬ 
bution to B(q = 0,cj). Our results obey the sum rule 
J 0 °° dujB(q,uj)(l + e _/3u; ) = D(q,r = 0). Additional sum 
rules are discussed in the appendix. 


IV. NUMERICAL RESULTS 

We first consider the spinful Holstein model [Eq. (1) 
with U = 0] with ujq = 0.5£, which exhibits a Peierls 
metal-insulator transition at A c ~ 0.25 [37, 38]. Figure 1 
shows the phonon spectral function B(q,uj) for different 


values of A. The system size and temperature were chosen 
as L = /3t = 42 (see also Sec. V). 

For a weak coupling A = 0.1, Fig. 1(a) reveals that the 
bare phonon mode at uj q = uj^ is renormalized near q = 0 
and <7 = 7 r, but largely unaffected at intermediate q. The 
small-g feature will be explained as a hybridization ef¬ 
fect in Sec. V. For a stronger coupling A = 0 . 2 , shown 
in Fig. 1(b), the renormalization of the phonon mode is 
significantly more pronounced, and we observe a partial 
softening near q = 7 r = 2 fcp, which is a precursor of the 
Peierls transition. Upon increasing A further, the phonon 
spectrum becomes gapless close to A c = 0.25. Beyond A c , 
the results in Figs. l(d)-l(f) are consistent (see below) 
with the emergence of a central (uj = 0) peak at q = 7 r, 
and a second peak at uj > 0. The central peak carries 
substantial spectral weight (which diverges with system 
size) and reflects the long-range lattice order with mod¬ 
ulation q = 2= 7 r. The phonon dispersion throughout 
the Brillouin zone appears to be most pronounced for 
A ~ A c , and becomes flatter in the Peierls phase. 

To understand the evolution of B(q, uj) at the ordering 
wavevector q = 7 r in more detail, we show in Fig. 2 the 
spectrum B(tt,uj) for the same parameters as in Fig. 1 . 
Except for A = 0.275, we observe two separate low-energy 
peaks in the phonon spectrum, although they are difficult 
to discern in the density plots of Fig. 1 . Starting at weak 
coupling A = 0.1, there is a peak at uj ~ 0 . 2 £ (correspond¬ 
ing to the finite-size charge gap in the dynamic charge 
structure factor at q = 7r, which scales as 1/L), and a 
peak close to the bare phonon frequency uj = ujq = 0.5£. 
With increasing A, both peaks move toward smaller uj. 
For the system size used, the spectrum becomes gapless 
close to the critical coupling A c « 0.25, coinciding with 
the appearance of the central peak. Figures 2(e) and 2 (f) 
reveal that a second peak emerges that becomes harder 
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FIG. 2. (Color online) Phonon spectral function B(q,u) at 
q = 7r for the same parameters as in Fig. 1. The dashed verti¬ 
cal lines indicate cjo = 0.5£. As in Fig. 1, we use a logarithmic 
scale and a plot range [0.1,10 3 ]. 


with increasing A but has small spectral weight compared 
to the central peak. 

It is also interesting to consider the impact of Coulomb 
repulsion in the framework of the Holstein-Hubbard 
model. In Fig. 3(a) we show the phonon spectral function 
for A = 0.3 and U = t, which should be compared to the 
results of Fig. 1(e). Quite remarkably, we find a com¬ 
parable renormalization of the bare phonon mode, but 
no central peak. The absence of the latter is expected 
since U = t is sufficient to drive the system from the 
Peierls into the intermediate metallic phase, cf. Fig. 8(a) 
in Ref. [37]. On the other hand, the renormalization 
of the phonon mode is much stronger than expected 
based on a simple effective Hubbard model with inter¬ 
action C/ e ff = U — 4£A, which is justified for sufficiently 
large uoo/t. For the parameters of Fig. 3(a), we have 
E/eff = —0.2 1. However, the phonon spectrum is renor¬ 
malized more strongly than in Fig. 1(a), where A = 0.1 
and E/ e ff = — 0.4t, suggesting that not only the effec¬ 
tive interaction but also the bare electron-phonon cou¬ 
pling determines the phonon renormalization. Finally, 
for U = 2t [Fig. 3(b)], corresponding to the Mott phase 
[37], the renormalization of the phonon mode is signifi¬ 
cantly weaker, and the phonon spectrum resembles quite 
closely the weak-coupling results shown in Fig. 1(a). 



FIG. 3. (Color online) Phonon spectral function B(q,uj) of 
the Holstein-Hubbard model. The dashed line corresponds to 
(jj 0 — 0.5 1. Here, A — 0.3, and /3t = L = 42. 

V. DISCUSSION 

In this section, we relate our numerical results to pre¬ 
vious work and to theoretical expectations. Additionally, 
we provide an explanation of the small-g feature observed 
in the metallic phase, and show how the discrete symme¬ 
try spontaneously broken at the Peierls transition can be 
captured in a continuum field theory. 

A. Relation to the charge structure factor 

As illustrated by the exact relation (7), the phonon 
Green function and the time-displaced charge correlation 
function are related by a convolution. Consequently, the 
phonon spectral function B(q, uj) and the dynamic charge 
structure factor 

N(q,u) = 4 \(m\p g \n)\ 2 e~ pEm 6(w - A nm ), (9) 

n,m 

with p q = L -1 / 2 ^2 r e zqr p ri contain in principle the same 
information, although the spectral weights may differ by 
orders of magnitude. In particular, as discussed in more 
detail in the Appendix, the dynamic charge structure 
factor also reveals the renormalized phonon excitations 
[21, 28, 31], and the phonon spectral function also con¬ 
tains signatures of the particle-hole continuum. 

Importantly, for values of q where the bare phonon 
mode and the particle-hole continuum overlap, the renor¬ 
malized phonon mode uj q will be damped, in contrast to 
the (^-function contribution suggested by the approxima¬ 
tion for B(q,w) given in Eq. (5) of Ref. [23]. For inter¬ 
mediate phonon frequencies ujq ~ £, no clear separation 
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exists between the phonon mode and particle-hole exci¬ 
tations. Finally, for ujq t (that is, larger than the bare 
bandwidth 4 1 of the particle-hole continuum), or deep in 
the Peierls phase where the lower edge of the particle- 
hole continuum lies above c5 g , we expect the phonon 
mode to remain a separate and well-defined excitation 
near q = 2/cp, as suggested by Figs. 3(b) and 2(a) in 
Ref. [23], respectively. 

An important corollary of the relation between B(q, c j) 
and N(q,w) concerns the shape of the spectrum in the 
vicinity of q = tt. In the Peierls phase, the size of the 
unit cell (Brillouin zone) doubles (halves). Although a 
perfect symmetry (extending to the spectral weights of 
excitations) between q = 0 and q = tt is only achieved in 
the limit A oo, this doubling implies that the renor¬ 
malized phonon frequency uj q cannot extrapolate to zero 
at q = 7r for A > A c . Such a gapless mode would neces¬ 
sarily have a counterpart near q = 0, and also in N(q, uS). 
However, gapless, long-wavelength particle-hole excita¬ 
tions are not compatible with an insulating Peierls state. 
We will see below that a simple low-energy theory in¬ 
stead suggests a gapped renormalized phonon mode and 
an isolated central peak at q = tt. This picture is consis¬ 
tent with our numerical results. On the other hand, we 
attribute the apparent existence of a gapless mode near 
q = 7r in results for the spinless Holstein model [21] to a 
failure to resolve the expected two-peak structure in the 
dynamic charge structure factor. 


B. Origin of the small-# feature 

The coupling between the bare phonon mode and the 
particle-hole excitations of 7 V(#,cj) [see Fig. 4(a)] via the 
electron-phonon interaction provides an explanation of 
the small-# feature visible in Fig. 1(a), as well as in the 
renormalized phonon frequency uj q in Refs. [17, 25]. In 
the analytical results of Refs. [17, 25], this feature occurs 
at a nonzero q and is quite sharp in momentum space, 
making it difficult to resolve fully in our finite-size data 
and explaining its absence in previous exact diagonaliza- 
tion results for small clusters [19]. 

Here, we explain this feature in terms of the hybridiza¬ 
tion between the bare phonon frequency uo q = ujq and 
particle-hole excitations that—in one dimension—have 
the linear dispersion uj = v^q for small #, see Fig. 4(a). 
The hybridization is captured by the random phase ap¬ 
proximation for the phonon propagator, 

D-\q, m n ) = [D°(*f2 n )] _1 + A 2 X °(q, *O n ), (10) 

where Q n is a bosonic Matsubara frequency, and 
X°(#, ifln) is the noninteracting charge susceptibility. 
Taking A as a free parameter to match the numerical re¬ 
sults, we find a hybridization at the intersection point of 
the free phonon dispersion and the free particle-hole ex¬ 
citations, as well as a gapless linear mode below the bare 
dispersion that corresponds to long-wavelength charge 
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FIG. 4. (Color online) (a) Dynamic charge structure factor of 
the noninteracting system (A — U = 0). (b) Phonon spectral 
function from Eq. (10) for A = 0.2. The dashed lines indicate 
ujq = 0.5t. Here, T — 0, and L — 402. 


fluctuations. While the CT-INT results in, for exam¬ 
ple, Fig. 1(a), do not fully resolve the hybridization, the 
agreement is satisfactory. Note that the random phase 
approximation does not take into account any phonon 
softening related to charge order, or the renormalization 
of v-p by interactions. The hybridization of charge and 
phonon modes at small q also follows from a Tomonaga- 
Luttinger model [46]. 

In accordance with the analytical results of Refs. [17, 
25], a hybridization of the phonon mode with the particle- 
hole continuum is also visible near q = 7r, and gives rise 
to damping of the phonon excitations. The above expla¬ 
nation suggests that the small-# feature is absent deep in 
the insulating Peierls phase, because the latter does not 
have low-energy excitations near q = 0. Accordingly, a 
suppression of the hybridization feature with increasing 
A is visible in Fig. 1. Finally, within the random phase 
approximation, the hybridization does not explain the 
observed hardening of the phonon mode near the zone 
boundary for ujq = 4 1 [19]. 


C. Finite-size effects 

Although our system size of L = 42 is significantly 
larger than in previous work, it is important to separate 
finite-size effects from generic features. 

The most notable finite-size artifact in our results is 
the charge gap at q = tt in the dynamic charge struc¬ 
ture factor and (for A > 0) also in the phonon spectral 
function. In the noninteracting case, this gap scales as 
1/L. While it is negligibly small in Fig. 4(a) {L = 402), 
it is about 0.3 1 for the system size L = 42 used in our 
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simulations. The results for the phonon spectral function 
in Fig. 2 reveal that with increasing A, the charge gap is 
reduced. At the critical point, the central peak at uj = 0 
appears. The closing of the finite-size charge gap almost 
exactly at the critical coupling A c is not expected to be 
generic, but to depend on system size and parameters. 

While a charge gap is expected in B(tt, uj) in the Peierls 
phase even in the thermodynamic limit, both B(q , uj) and 
N(q,ui) are gapless at q = tt throughout the metallic 
phase for L oo [23]. In the insulating Peierls phase, 
the central peak at u = 0 is separated from excitations 
with uj > 0 by the interaction-generated charge gap. Im¬ 
portantly, the lower edge of the particle-hole continuum 
in the metallic phase corresponds to a branch cut (simi¬ 
lar to other excitations of ID systems) with finite spec¬ 
tral weight, whereas the weight of the central peak at 
uo = 0 in the Peierls phase diverges in the thermody¬ 
namic limit. Finite-size effects on the phonon spectrum 
of a spin-Peierls model have been discussed in Ref. [47]. 

To avoid spurious finite-size effects, the system size 
should be large enough to have a noninteracting charge 
gap smaller than the bare phonon frequency. Other¬ 
wise, there is no coupling between the bare mode and 
the particle-hole continuum, which is not generic for the 
adiabatic regime. In the latter, and more generally for 
any c < 4£, the bare phonon mode lies on top of the 
noninteracting particle-hole continuum in extended re¬ 
gions of the Brillouin zone. A nonzero electron-phonon 
coupling then gives rise to a mixing of these excitations, 
leading to renormalization and a finite lifetime of phonon 
excitations. In this sense, Eq. (5) in Ref. [23] should be 
regarded as an approximate result because it suggests 
the existence of phonon excitations with infinite lifetime 
described by S(cj — ui q ). 


D. Soft-mode versus central-peak behavior 

Peierls transitions are often classified as either soft- 
mode or central-peak transitions. In terms of the cou¬ 
pling A (rather than temperature, as appropriate for ex¬ 
periments on quasi-ID systems), a ID soft-mode transi¬ 
tion involving a single phonon mode is characterized by 
a softening uo 2 k F —> 0 for A < A c , a completely soft mode 
uj 2 k F = 0 at A = A c , and a subsequent hardening (i.e., 
co> 2 k F > 0) in the Peierls phase with an additional central 
peak at u = 0 reflecting the long-range lattice order. In 
contrast, for a central-peak transition, 0J2k F stays nonzero 
or even hardens across the transition, and a central peak 
appears at the critical point and persists for A > A c . 

The above simple picture is modified in several ways. 

(i) Because of the relation to the charge structure factor, 
B(q,uj) in general has a continuum of excitations. Al¬ 
though the spectra shown here for the adiabatic regime 
are dominated by a few peaks, the particle-hole contin¬ 
uum has substantial spectral weight for ujo > t [19, 23]. 

(ii) In the thermodynamic limit, B(q,uj) has soft (i.e., 
gapless) excitations in the metallic phase corresponding 


to the gapless particle-hole continuum at q = tt. (iii) In 
the Peierls phase, the renormalized phonon mode (which 
has a finite energy for A > A c ) has very small spectral 
weight compared to the divergent central peak, which 
may make its identification difficult for both numerical 
simulations and experiments. Similar behavior has been 
reported for a spin-Peierls model in Ref. [47]; in particu¬ 
lar, it was pointed out that the two peaks in B(q, tt) are 
only visible on large enough system sizes. Finally, it is 
not clear if the hardening of uj 2 k F in the Peierls phase can 
be clearly distinguished from the opening of the gap in 
the particle-hole continuum in the vicinity of A c . 

Keeping in mind these complications, our numerical 
results in Figs. 1 and 2 are nevertheless compatible with 
a soft mode transition. The dispersive feature that ex¬ 
ists throughout the Brillouin zone and emerges from the 
bare phonon mode uj q = Co>o upon turning on the electron- 
phonon interaction can be identified with the renormal¬ 
ized phonon frequency uj q . We observe that Lo q soft¬ 
ens near q = 7r as A —)> A c , and hardens again in the 
Peierls phase. However, due to limitations in system size 
and spectral resolution, we are unable to unambiguously 
demonstrate the existence of a true soft mode at A c . 

Let us briefly discuss relevant previous works on this 
issue. In Refs. [17, 23, 25], the renormalized phonon 
frequency uj q of the spinless Holstein model was found 
to soften at the critical point, and harden again in the 
Peierls phase. This behavior is consistent with a soft- 
mode transition, but the results for the phonon spectral 
function in Fig. 2 of Ref. [23] do not seem to show a 
central peak related to long-range order, in contrast to 
numerical results [19]. A softening is also visible in ex¬ 
act diagonalization results for small clusters in Ref. [19]. 
Based on QMC simulations and fits of the phonon prop¬ 
agator to a simple two-mode form, a softening (harden¬ 
ing) in the metallic (insulating) phase was suggested in 
Ref. [22]. Although we can reproduce these results, it 
is not clear how reliable such fits are given the large ra¬ 
tio of spectral weights between the central peak and the 
phonon peak. Moreover, we find that the minimum of the 
fitted phonon frequency does not give the correct critical 
value A c in the adiabatic regime (cjo = 0.1 1). Finally, 
the phonon softening and ensuing hardening observed in 
Migdal-Eliashberg theory [48] and in approximate solu¬ 
tions of the Holstein model [49] for any band filling ap¬ 
pear to be a consequence of the breakdown of the adia¬ 
batic approximation. A soft phonon mode seems to be 
always linked with a lattice instability (phase transition), 
which requires a commensurate band filling. 


E. Field-theory description 

The observation of a Peierls transition implies the ex¬ 
istence of a Z 2 symmetry which, in one dimension, can 
be spontaneously broken at zero temperature and would 
account for the observed central peak. Here we provide a 
minimal continuum theory that reveals such a Z 2 symme- 
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try for the Holstein model. For simplicity, we consider 
spinless fermions and assume an arbitrary value of the 
Fermi wave vector kp. 

Near dzfcp, the electronic degrees of freedom can be 
described by a linear dispersion, leading to the Dirac form 

H e = V F J2^-*t<Tz* k , ( 11 ) 

k 


where \E^ = is a spinor of creation operators 

for left and right moving fermions, respectively, and cr z 
is the usual Pauli matrix. H e has a Ug z (1) symmetry 
related to the pseudospin S z = \ J2k we have 

Ue^kUe = U e ^ k with U e = and U e = 

The free phonon part can be written as 

k 

with ujk = co>o for optical phonons. 

In the framework of bosonization, the Peierls transition 
arises from phonon-mediated single-particle backscatter- 
ing described by (see Appendix B of Ref. [31]) 

Hi = ^2, \n\L k _q (o^_2k F - q + ^2/cp+g) 

+ Lk-qRk (^2/c F +g + ®-2k F -q) ' 

To derive a low-energy theory, we only consider small 
momenta \q\ <C 2 kp. Consequently, for our purposes, 
<^- 2 k F -q an( l ^2/cp+g ma y ^e re g ar( ied as independent 

bosonic modes [50]. With a J = (d^ 2/cp _ g , ^\k F +q) Ihe 
phonon Hamiltonian can be written as 

H p = ^2at\[u + (q)+u-(q)a z \a q , (14) 

q 

where u±(q) = (uj- 2 k F - q ± uo 2 k F +q) /2. Under the above 
assumption of two independent modes, H p has two 17(1) 
symmetries generated by the conservation of the total 
phonon number N = J2 q a^ q a q , and of N z = ^2 q a^ q a z a q ^ 
denoted as Un(1) and Un z (1), respectively. 

The electron-phonon interaction term (13) breaks the 
full 17 jv( 1) x Un z { 1) x Us z ( 1) of the free Hamilto¬ 
nian (in particular, the number of phonons is not con¬ 
served). However, under a combined transformation 
U^o = e l( t >N z e ies z the interaction term transforms as 


u; 


:,lHi % <e = -)= Y^\uta-U g ^ k _ q e^ (15) 

V k,q 

X (^-2fe P -q + “2fe F + ? ) + H.C. , 


where a± = \ ( a x =b icr y ). Since UqCT-Uq = e 2( V_, set¬ 
ting 0 = —(f) shows that S z — N z is a generator of a 
Us z -n z { 1) symmetry of the full Hamiltonian H e + H p + 
H\. On the other hand the interaction term reduces the 


Un{ 1) symmetry to a discrete Z 2 symmetry. Therefore, 
the symmetry of the full Hamiltonian is Us z -n z { 1) x Z 2 . 

As mentioned previously, the Us z ~n z (1) x Z 2 sym¬ 
metry can explain the existence of a Peierls state with 
long-range order at T = 0. The Us z -n z { 1) symmetry ac¬ 
counts for a combined charge and phonon mode, whereas 
the Z 2 symmetry, if broken at T = 0, gives rise to the 
central peak feature. Within our continuum field-theory 
description, commensurate band fillings cannot be un¬ 
ambiguously defined. Nevertheless, the reduction of the 
C/jv(1) symmetry to Z 2 accounts for the expected pin¬ 
ning of the charge-density wave for commensurate band 
fillings [51]. The Us z -n z { 1) symmetry is reminiscent of 
the C7(l) Gross-Neveu theory with one bosonic mode and 
two fermion flavors [52] that describes, for example, the 
transition from a Dirac semimetal to a Kekule ordered 
state in two dimensions [53]. 


VI. CONCLUSIONS 

We showed that the phonon spectral function of 
electron-phonon models can be calculated using the CT- 
INT QMC method with the help of a generating func¬ 
tional. The same idea can be used to measure other cor¬ 
relation functions and observables, and for other models, 
which further extends the usefulness and versatility of 
the CT-INT method. 

Our results for the ID Holstein-Hubbard model are 
consistent with a soft-mode Peierls transition in the adi¬ 
abatic regime, characterized by a softening of the q = 2 kp 
phonon excitations, a gapless mode at the critical point, 
and a subsequent hardening of the phonon mode in the 
Peierls phase. However, this simple picture is compli¬ 
cated by finite-size effects, and by the mixing of phonon 
and charge excitations mediated by the electron-phonon 
interaction. We explained a small-g anomaly of the 
metallic phase, previously observed for the spinless Hol¬ 
stein model, in terms of a hybridization of the bare 
phonon mode and the particle-hole continuum. Overall, 
our results for the spinful Holstein model are very similar 
to previous results for the spinless Holstein model. 

For a Hubbard repulsion large enough to yield a metal¬ 
lic state, we observed a suppression of the central peak 
related to long-range order, but a significant remaining 
renormalization of the phonon mode that only disappears 
in the Mott phase for even larger Hubbard interaction. 

Finally, we provided a unified picture of current and 
previous work, and showed how the discrete symmetry 
spontaneously broken in the Peierls phase can be cap¬ 
tured in a continuum field theory. 
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Appendix: Exact relation between phonon and 
charge spectra 


gives an additional contribution at uj = ±co>o. Here, V 
denotes the principal value. 

To derive Eq. (A.4), we used a partial-fraction decom¬ 
position, leading to poles of both first and second or¬ 
der. However, poles of second order are forbidden by the 
Lehmann representation (A.l). Therefore, their weights 
have to vanish, which (for uj$ > 0) leads to the sum rule 


For the Holstein model, the phonon spectral function 
B(q,uj) and the dynamic charge structure factor N(q,uj) 
in principle contain the same information. Here, starting 
from Eq. (7), we derive an exact relation between the 
spectral functions as well as additional sum rules, and 
discuss the implications. 


1. Analytic properties of the spectral functions 

To simplify the notation, we define B(q,uj) = (1 — 
e~ /3uj )B(q, uj) and N(q,uj) = (1 — e~^ UJ )N(q, uj). The 
phonon spectral function B(q,uj) can be obtained from 
the Lehmann representation of the phonon propagator 

D{q,z) = - f duj (A.l) 

J — OO Z ^ 

by analyzing its pole structure in the complex-frequency 
plane. For simplicity, we restrict our considerations to 
finite Hilbert spaces, where D(q , z) has only simple poles 
on the real axis [54] determined by the exact relation 

D(q, z) = D°(z) + g 2 D°(z) 2 x(q, z). (A.2) 

The term ~ g 2 in Eq. (A.2) gives rise to a product of 
poles arising from the free phonon propagator D°(z) = 
— z 2 )- 1 and the charge susceptibility 

x(q,z) = (p q (z)p- q ) = - f du> . (A.3) 

J -oo Z-U 

A partial-fraction decomposition and comparison of the 
pole structure of the two sides of Eq. (A.2) gives 

B(q,u>) = B°(q,uj) + g 2 D 0 (u>) 2 N(q,u>) + B 1 (q,u). 

(A.4) 


For <7 = 0 , the phonon spectral function is given by 

B°(q, w) = ^T__ [£( w - w 0 ) - 6(u + w 0 )] , (A.5) 

which describes excitations at the bare phonon frequency 
u j = ±co>o. Any finite electron-phonon coupling leads to 
the appearance of two additional terms in Eq. (A.4): 
The first contains the whole charge spectrum N(q,uj) 
reweighted by the free phonon propagator, while 


B 1 (q,co) = [<5(w - w 0 ) - <5(w + w 0 )] 

ido 


,oo (A.6) 

xP dio 1 lj 1 D° (u>') 2 N (q, u>') 

Jo 


V du 2 W 2 N(q,u) =0. (A.7) 

Jo W ^0 

From Eqs. (A.4) and (A.7), we obtain an equivalent sum 
rule for 5(g,a;), 

duo uj (c S 2 — cJq) 5(g, uj) = 0 , (A. 8 ) 

which is just a combination of the first and third moment 
of B(<q,ui). In the same way, the absence of higher-order 
poles requires 7V(g,cj = Tcjo) = 0. 

For finite electron-phonon coupling, Eq. (A.2) can also 
be used to obtain the charge spectrum 

N(q,u) = ^ (u 2 -ul) 2 B(q,ui) (A.9) 

from the phonon spectral function. Here, contributions 
at uj = ±co>o are removed from N(q,uj) by the prefactor. 


2. Implications for the spectral properties 

According to Eqs. (A.4) and (A.9), H(g, uj) and N(q, uj) 
share the same spectral information, up to an additional 
contribution to B(q , ±co>o) that consists of the free phonon 
spectrum B°(q, uj) and a compensating term H x (g, uj) due 
to finite interactions. 

For q = 0, because of charge conservation, N(q, uj) only 
has a static contribution at uj = 0. Thus, B x (q = 0,cj) = 
0 and the full phonon spectrum is given by the free part 
at uj = ±co>o and the static contribution to N(q,uj). 

For q 7 ^ 0, any finite electron-phonon coupling seems to 
shift the phonon dispersion away from uj = Tcjo- Exact 
diagonalization data for the spinless Holstein model [19] 
suggest that B(q,±uj 0 ) vanishes and therefore B°(q,uj) 
and B 1 (q 1 uj) compensate each other [55]. In general, for 
q 7 ^ 0, both H(g, uj) and N(q, uj) contain signatures of the 
phonon dispersion as well as the particle-hole continuum, 
although the spectral weights may be very different. 

The condition B(q,uJo) > 0 sets an upper bound to the 
integral in Eq. (A. 6 ), 

roo 

V / duJ 7~2 -< 2-2 . (A.10) 

Jo (uj — uJq) z 9 

For the integral to converge, N(q,uj) has to vanish when 
approaching ujq. Thus, a nonzero electron-phonon inter¬ 
action splits the charge spectrum at uj = ujq. 

Further insight into the distribution of spectral weight 
can be obtained from the sum rules (A.7) and (A. 8 ). We 
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restrict our discussion to B(q,cj), but the same argu¬ 
ments hold for N(q,uo). For uo > 0, B(q,uo) > 0 but 
the prefactor ( [u 2 — uJq) changes sign at c^o- This sign 
change divides the frequency axis into regions uj < uoo and 
uj > cjo, whose integrated spectral weights have to com¬ 
pensate each other in the sum rule [56]. Note that spec¬ 
tral weight at uo = 0 and uj = ujq does not contribute to 
the sum rule, therefore the noninteracting phonon disper¬ 
sion fulfills Eq. (A.8) trivially. By adiabatically switching 
on the electron-phonon coupling, the particle-hole contin¬ 
uum enters B(q,uj) and spectral weight has to be redis¬ 


tributed to fulfill Eq. (A.8). For wave vectors such that 
the particle-hole continuum only enters one of the two re¬ 
gions, spectral weight has to appear in the other region. 
This can be most easily achieved by shifting the phonon 
dispersion. Both the hardening of the phonon disper¬ 
sion for cjo t [19, 23], and the hybridization with the 
particle-hole continuum as well as the phonon softening 
for cjo £, are consistent with the sum rule (A.8). Fur¬ 
thermore, in the Peierls phase, the charge gap (the lowest 
excitation at q = n) cannot become larger than ujq, as the 
central peak does not contribute to the sum rule (A.8). 
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